subroutine CalcPointsFacesCellsNeigForVisc()
use Data_Mesh
implicit none

    integer :: i,j, iPoint,jPoint, iFace,iCell, iCellPoint, iCellFace, jCell, iFacePoint, jFace
    integer :: nPoint2NeigCellsMax, nFace2NeigCellsMax, nCell2NeigCellsMax
    integer :: n,m, k, iCellFaceNeighbours(100), iCellNeigh, jCellFace, jFacePoint
    integer :: iCellOwner, iCellNeighbour
    logical :: newCell, nonewFace, NeighFace
    integer,allocatable :: nCellsForPoint(:)
    integer,allocatable :: nCellsForFace(:)
    integer,allocatable :: Face2NeigCellsTemp(:,:)
    integer,allocatable :: Cell2NeigCellsTemp(:,:)

    ! compute nPoint2NeigCells and Point2NeigCells
    allocate(nCellsForPoint(nPoints))
    
    nCellsForPoint = 0
    
    do iCell = 1,nCells
        do iCellPoint = 1,8
            iPoint = meshCellPoints(iCellPoint,iCell)
            nCellsForPoint(iPoint) = nCellsForPoint(iPoint) + 1
        end do
    end do
    nPoint2NeigCellsMax = maxVal(nCellsForPoint)
    deallocate(nCellsForPoint)
    
    allocate(nPoint2NeigCells(nPoints))
    allocate(Point2NeigCells(nPoint2NeigCellsMax,nPoints))
    
    nPoint2NeigCells = 0
    Point2NeigCells = 0
    do iCell = 1,nCells
        do iCellPoint = 1,8
            iPoint = meshCellPoints(iCellPoint,iCell)
            newCell = .true.
            do i = 1,nPoint2NeigCells(iPoint)
                if(Point2NeigCells(i,iPoint).eq.iCell)then
                    newCell = .false.
                    exit
                end if
            end do
            if(newCell)then
                nPoint2NeigCells(iPoint) = nPoint2NeigCells(iPoint) + 1
                Point2NeigCells(nPoint2NeigCells(iPoint),iPoint) = iCell
            end if
        end do
    end do

    ! ************************************************** !

    ! compute nFace2NeigCells and Face2NeigCells
    allocate(Face2NeigCellsTemp(100,nFaces))
    allocate(nFace2NeigCells(nFaces))

    Face2NeigCellsTemp = 0
    nFace2NeigCells = 0

    do iFace = 1,nFaces
        do iFacePoint = 1,4
            iPoint = meshFaces(iFacePoint,iFace)
            do i = 1,nPoint2NeigCells(iPoint)
                jCell = Point2NeigCells(i,iPoint)
                newCell = .true.
                do j = 1,nFace2NeigCells(iFace)
                    if(Face2NeigCellsTemp(j,iFace).eq.jCell)then
                        newCell = .false.
                        exit
                    end if
                end do
                if(newCell)then
                    nFace2NeigCells(iFace) = nFace2NeigCells(iFace) + 1
                    Face2NeigCellsTemp(nFace2NeigCells(iFace),iFace) = jCell
                end if
            end do
        end do
    end do

    nFace2NeigCellsMax = maxVal(nFace2NeigCells)
    allocate(Face2NeigCells(nFace2NeigCellsMax,nFaces))
    Face2NeigCells = 0
    
    do iFace = 1,nFaces
        do i = 1,nFace2NeigCells(iFace)
            Face2NeigCells(i,iFace) = Face2NeigCellsTemp(i,iFace)
        end do
    end do

    deallocate(Face2NeigCellsTemp)
    ! ************************************************** !
!    allocate(nCellsForFace(nFaces))
!    
!    nCellsForFace = 0
!    
!    do iCell = 1,nCells
!        do iCellFace = 1,6
!            iFace = meshCellPoints(iCellFace,iCell)
!            nCellsForFace(iFace) = nCellsForFace(iFace) + 1
!        end do
!    end do
!    nFace2NeigCellsMax = maxVal(nCellsForFace)
!    deallocate(nCellsForFace)
!
!    allocate(nFace2NeigCells(nFaces))
!    allocate(Face2NeigCells(nFace2NeigCellsMax,nFaces))
!    
!    nFace2NeigCells = 0
!    Face2NeigCells = 0
!    do iCell = 1,nCells
!        do iCellFace = 1,6
!            iFace = meshCellPoints(iCellFace,iCell)
!            newCell = .true.
!            do i = 1,nFace2NeigCells(iFace)
!                if(Face2NeigCells(i,iFace).eq.iCell)then
!                    newCell = .false.
!                    exit
!                end if
!            end do
!            if(newCell)then
!                nFace2NeigCells(iFace) = nFace2NeigCells(iFace) + 1
!                Face2NeigCells(nFace2NeigCells(iFace),iFace) = iCell
!            end if
!        end do
!    end do
    ! ************************************************** !

    ! compute nCell2NeigCells and Cell2NeigCells
    allocate(Cell2NeigCellsTemp(100,nCells))
    allocate(nCell2NeigCells(nCells))

    Cell2NeigCellsTemp = 0
    nCell2NeigCells = 0

    do iCell = 1,nCells
        do iCellPoint = 1,8
            iPoint = meshCellPoints(iCellPoint,iCell)
            do i = 1,nPoint2NeigCells(iPoint)
                jCell = Point2NeigCells(i,iPoint)
                newCell = .true.
                do j = 1,nCell2NeigCells(iCell)
                    if(Cell2NeigCellsTemp(j,iCell).eq.jCell)then
                        newCell = .false.
                        exit
                    end if
                end do
                if(newCell)then
                    nCell2NeigCells(iCell) = nCell2NeigCells(iCell) + 1
                    Cell2NeigCellsTemp(nCell2NeigCells(iCell),iCell) = jCell
                end if
            end do
        end do
    end do

    nCell2NeigCellsMax = maxVal(nCell2NeigCells)
    allocate(Cell2NeigCells(nCell2NeigCellsMax,nCells))
    Cell2NeigCells = 0
    
    do iCell = 1,nCells
        do i = 1,nCell2NeigCells(iCell)
            Cell2NeigCells(i,iCell) = Cell2NeigCellsTemp(i,iCell)
        end do
    end do

    deallocate(Cell2NeigCellsTemp)
    ! ************************************************** !

! ***************************************************************************************** !
! ***************************************************************************************** !

!    allocate(CellStencil(nCells))
!
!    do iCell = 1,nCells
!
!        CellStencil(iCell)%nCells = nCell2NeigCells(iCell)
!        n = CellStencil(iCell)%nCells
!
!        allocate(CellStencil(iCell)%iCell(n))
!
!        do i = 1,n
!            CellStencil(iCell)%iCell(i) = Cell2NeigCells(i,iCell)
!        end do
!    end do
!
!! ****************************************************** !
!
!    do iCell = 1,nCells
!        CellStencil(iCell)%nFaces = 6
!
!        do i = 1,6
!            iCellFaceNeighbours(i) = meshCells(i,iCell)
!        end do
!
!        n = CellStencil(iCell)%nCells
!        do iCellNeigh = 1,n
!            jCell = CellStencil(iCell)%iCell(iCellNeigh)
!            if(jCell.eq.iCell) cycle
!            do jCellFace = 1,6
!                jFace = meshCells(jCellFace,jCell)
!                
!                noNewFace = .false.
!                do i = 1,CellStencil(iCell)%nFaces
!                    if(iCellFaceNeighbours(i).eq.jFace)then
!                        noNewFace = .true.
!                        exit
!                    end if
!                end do
!
!                if(noNewFace) cycle
!
!                NeighFace = .false.
!                do jFacePoint = 1,4
!                    jPoint = meshFaces(jFacePoint,jFace)
!                    do iCellPoint = 1,8
!                        iPoint = meshCellPoints(iCellPoint,iCell)
!                        if(iPoint.eq.jPoint)then
!                            NeighFace = .true.
!                            exit
!                        end if
!                    end do
!                end do
!                
!                if(NeighFace)then
!                    CellStencil(iCell)%nFaces = CellStencil(iCell)%nFaces + 1
!                    i = CellStencil(iCell)%nFaces
!                    iCellFaceNeighbours(i) = jFace
!                end if
!            end do
!        end do
!        
!        n = CellStencil(iCell)%nFaces
!        allocate(CellStencil(iCell)%iFace(n))
!        do i = 1,n
!            CellStencil(iCell)%iFace(i) = iCellFaceNeighbours(i)
!        end do
!        
!        ! define owner and neighbour cell
!        allocate(CellStencil(iCell)%FaceOwner(n))
!        allocate(CellStencil(iCell)%FaceNeighbour(n))
!        allocate(CellStencil(iCell)%InternalFace(n))
!        CellStencil(iCell)%FaceOwner = 0
!        CellStencil(iCell)%FaceNeighbour = 0
!        CellStencil(iCell)%InternalFace = .false.
!
!        do i = 1,n
!            iFace = CellStencil(iCell)%iFace(i)
!            iCellOwner = meshFaceOwner(iFace)
!            if(iFace.gt.nInternalFaces)then
!                CellStencil(iCell)%InternalFace(i) = .true.
!                iCellNeighbour = 0
!            else
!                iCellNeighbour = meshFaceNeighbour(iFace)
!            end if
!            
!            do j = 1,CellStencil(iCell)%nCells
!                jCell = CellStencil(iCell)%iCell(j)
!                if(iCellOwner.eq.jCell)then
!                    CellStencil(iCell)%FaceOwner(i) = j
!                end if
!                if(iCellNeighbour.eq.jCell)then
!                    CellStencil(iCell)%FaceNeighbour(i) = j
!                end if
!            end do
!        end do
!
!    end do
!
!! nFace2Point  Face2Point
!    do iCell = 1,nCells
!        n = CellStencil(iCell)%nFaces
!        allocate(CellStencil(iCell)%nFace2Point(n))
!        allocate(CellStencil(iCell)%Face2Point(8,n))
!
!        CellStencil(iCell)%nFace2Point = 0
!        CellStencil(iCell)%Face2Point  = 0
!
!        do i = 1,n
!            iFace = CellStencil(iCell)%iFace(i)
!            do j = 1,4
!                jPoint = meshFaces(j,iFace)
!                do k = 1,8
!                    iPoint = meshCellPoints(k,iCell)
!                    if(iPoint.eq.jPoint)then
!                        CellStencil(iCell)%nFace2Point(i) = CellStencil(iCell)%nFace2Point(i) + 1
!                        m = CellStencil(iCell)%nFace2Point(i)
!                        CellStencil(iCell)%Face2Point(m,i) = k
!                    end if
!                end do
!            end do
!        end do
!    end do

end subroutine CalcPointsFacesCellsNeigForVisc
